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ABSTRACT 


Massive disk galaxies like our Milky Way should host an ancient, metal-poor, and centrally concen- 
trated stellar population. This population reflects the star formation and enrichment in the few most 
massive progenitor components that coalesced at high redshift to form the proto-Galary. While metal- 
poor stars are known to reside in the inner few kiloparsecs of our Galaxy, current data do not yet provide 
a comprehensive picture of such a metal-poor “heart” of the Milky Way. We use information from 
Gaia DR3, especially the XP spectra, to construct a sample of 2 million bright (Ggp < 15.5 mag) giant 
stars within 30° of the Galactic Center with robust [M/H] estimates, 6[M/H] < 0.1. For most sample 
members we can calculate orbits based on Gaia RVS velocities and astrometry. This sample reveals an 
extensive, ancient, and metal-poor population that includes ~ 18, 000 stars with —2.7 < [M/H] < —1.5, 
representing a stellar mass of > 5 x 10’Mo. The spatial distribution of these [M/H] < —1.5 stars 
has a Gaussian extent of only org, ~ 2.7 kpc around the Galactic center, with most of these or- 
bits being confined to the inner Galaxy. At high orbital eccentricities, there is clear evidence for 
accreted halo stars in their pericentral orbit phase. Stars with [M/H] < —2 show no net rotation, 
whereas those with [M/H] ~ —1 are rotation dominated. Most of the tightly bound stars show [a/Fe]- 
enhancement and [Al/Fe|-[Mn/Fe] abundance patterns expected for an origin in the more massive 
portions of the proto-Galaxy. These central, metal-poor stars most likely predate the oldest part of 
the disk (Tage © 12.5 Gyrs), which implies that they formed at z 2 5, forging the proto-Milky Way. 
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1. INTRODUCTION 


Understanding the formation history of our own 
Galaxy, especially its earliest phases, has been a cen- 
tral goal of Galactic Archeology for decades (Freeman 
& Bland-Hawthorn 2002), with dramatic progress en- 
abled by a suite of spectroscopic surveys (e.g., Yanny 
et al. 2009; De Silva et al. 2015; Majewski et al. 2017; 
Conroy et al. 2019) and ESA’s Gaia mission (Perryman 
et al. 2001; Gaia Collaboration et al. 2022). The earliest 
phases of the Milky Way’s star-formation and enrich- 
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ment history are reflected at the present epoch in the 
orbit- and abundance- distribution of old and metal- 
poor stars. In the context of the hierarchical forma- 
tion of massive disk galaxies like the Milky Way, we 
should expect the oldest and most metal-poor stars (say 
[M/H] < —1.5) to be a mix of stars that a) formed within 
one of the main overdensities which coalesced early to 
form the proto-Galaxy, or b) formed early in distinct 
satellite galaxies that eventually merged with the main 
body (Zolotov et al. 2009). The first channel is com- 
monly referred to as in situ formation, the second ac- 
creted (Tumlinson 2010; Pillepich et al. 2015; El-Badry 
et al. 2018; Renaud et al. 2021b). If the in situ stars 
formed in a considerably more massive potential well 
than the accreted stars, this difference in origin should 
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also be reflected in their stars’ abundance patterns, such 
as [a/Fe] vs. [Fe.H] or [Mg/Mn] vs. [Al/Fe] (Zolotov 
et al. 2010; Hawkins et al. 2015; Das et al. 2020; Horta 
et al. 2021; Belokurov & Kravtsov 2022; Conroy et al. 
2022, and others). 

However, at very early epochs, this distinction be- 
tween in situ and accreted may become rather blurry. 
High-resolution formation simulations of galaxies like 
the Milky Way show (e.g. Renaud et al. 2021a,b) that 
pieces of comparable mass may rapidly coalesce early on 
(say, z > 4) in a sequence of “major mergers”. One ter- 
minology choice is to label only the one piece that was 
somewhat more massive as the in situ component, and 
to label all other pieces as accreted (which put together 
may constitute the majority of mass). However, simula- 
tions show that it may be operationally impossible later 
on to differentiate what was in situ and what accreted at 
very early epochs, whether one uses orbit or abundance 
information (e.g. Renaud et al. 2021b; Brauer et al. 2022; 
Orkney et al. 2022). Alternatively, one could label col- 
lectively all the major pieces that coalesced very early 
on, say at z 2 5, as the proto-Galazy. (see Conroy et al. 
2022)', and then differentiate subsequent additions at 
slightly later epochs (say, z < 3) as either in situ, if 
the material was brought in as gas, or as accreted, if 
the stars formed in a distinct potential well that then 
subsequently merged. Here, we opt for the second ter- 
minology, in part in light of the results we find. Conse- 
quently, we refer to the oldest parts of the Milky Way 
as the proto-Galazy, without trying to single out one of 
the contributing pieces as in situ. 

The nature of the accreted component of our Galaxy 
at Rec 2 5 kpc has come into focus due to the com- 
bination of Gaia and large ground-based spectroscopic 
surveys. Tidal debris from the disrupted Sagittarius 
satellite (Ibata et al. 1994) dominates the halo from 
20 — 50 kpc with stars mostly on highly inclined orbits 
with substantial angular momentum (Majewski et al. 
2003; Naidu et al. 2020); at 5 — 25 kpc the debris from 
the disrupted Gaia-Sausage-Enceladus (GSE) satellite 
(Helmi et al. 2018; Belokurov et al. 2018b; Naidu et al. 
2020) dominates, with stars on highly eccentric orbits 
(e = 0.7) with little (slightly retrograde) angular mo- 
mentum. Beyond these dominant components, a grow- 
ing number of additional, distinct accreted components 
of the Galaxy have been identified (e.g., Newberg et al. 
2009; Kruijssen et al. 2019; Myeong et al. 2019; Naidu 
et al. 2020; Yuan et al. 2020; Malhan et al. 2022). Be- 
yond their orbits, the abundance patterns of these stars 
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also point toward an accreted origin: their distribu- 
tion in the [a/M]-[M/H] plane (lower [a/M] at a given 
[M/H]) indicates that their birth material was enriched 
in a potential well of lower mass, i.e. a satellite galaxy 
(e.g. Frebel & Norris 2015; Hawkins et al. 2015; Lee et al. 
2015). In addition to presumed satellite galaxy debris, 
the stellar halo also appears to encompass a large num- 
ber of disrupted star clusters (e.g. Malhan et al. 2018; 
Shipp et al. 2018; Bonaca et al. 2021). 

There has also been progress in understanding the old 
proto-Galactic or in situ component of the Milky Way. 
The old and metal-poor Galactic in situ component that 
has been mapped best is the old, a-enhanced disk (e.g., 
Hayden et al. 2015; Bonaca et al. 2020; Belokurov et al. 
2020; Xiang & Rix 2022): its stars seem to date back to 
> 12.5 Gyrs (or z 2 5), are centrally concentrated and 
form a thick disk. However, the old a-enhanced disk 
metallicity distribution function (MDF) is truncated, or 
at least drops sharply, below [M/H] = —1. It seems in- 
evitable that there must have been a substantive proto- 
Galactic stellar population responsible for enriching the 
old disk’s birth material to [M/H] = —1. Recent surveys 
have provided clear evidence for a so-called in situ halo 
component in the Milky Way (e.g., Bonaca et al. 2017; 
Haywood et al. 2018; Di Matteo et al. 2019; Naidu et al. 
2020; Bonaca et al. 2020; Belokurov et al. 2020; Horta 
et al. 2021). Many of these studies drew on samples 
at the Solar radius and beyond, and those in situ halo 
stars may reflect very early disk stars that were kicked 
up, or splashed to hotter orbits. 

While the chemodynamical structure of the Galaxy 
at Rac 2 5 kpc has come into focus over the past 
decade, the nature of the inner Galaxy — especially at 
low metallicity — has proven more elusive. Fortunately, 
recent work has begun to piece together a picture of the 
metal-poor inner Galaxy with limited sample sizes (e.g., 
Ness et al. 2013; Garcia Pérez et al. 2013; Schlaufman 
& Casey 2014; Casey & Schlaufman 2015; Ness et al. 
2015; Koch et al. 2016; Garcia Pérez et al. 2018; Reg- 
giani et al. 2020; Arentsen et al. 2020a; Lucey et al. 
2021). In particular, narrow-band photometric surveys 
have spectacularly enabled a more efficient selection of 
metal-poor stars towards the inner Galaxy (e.g., Ar- 
entsen et al. 2020a,b), providing metal-poor samples in 
the inner Galaxy of several thousand objects. Spectro- 
scopic follow-up has shown that these stars are indeed 
metal poor and show very little rotation at the lowest 
metallicities (Arentsen et al. 2020b). Kruijssen et al. 
(2019, 2020) and Forbes (2020) recently used the proper- 
ties of inner Galaxy globular clusters as tracers of a very 
old Milky Way component, which they dubbed Kraken 
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and Koala. They considered these possible instances of 
early accretion. 

Most recently, Belokurov & Kravtsov (2022) and Con- 
roy et al. (2022) have been able to shed first light on 
the question of how the old disk emerged from an ear- 
lier, more metal-poor population of stars on dynami- 
cally hot orbits. Belokurov & Kravtsov (2022) combined 
APOGEE and Gaia astrometry to study the kinemat- 
ics of stars in the inner Galaxy down to metallicities of 
[M/H] ~ —1.5. They found that metal-poor stars in the 
inner Galaxy with abundance patterns they identified 
with in situ formation (dubbed Aurora) show rather lit- 
tle net rotation, and reasoned that these stars preceded 
the Milky Way’s ancient disk (see also Arentsen et al. 
2020a). Conroy et al. (2022) used chemistry, kinemat- 
ics, and ages from H3 and Gaia data to identify likely 
proto-Galactic stars down to [M/H] ~ —2.5 with ages 
2 13 Gyr. This recent work clearly implies that our 
Milky Way has some form of a proto-Galactic popula- 
tion. Whether or not there exists an operationally sep- 
arable accreted population buried deep in the potential 
well of our Galaxy remains an open question (cf. Krui- 
jssen et al. 2020; Horta et al. 2021; Myeong et al. 2022). 

It is on this background that we aim to flesh out 
a more comprehensive picture of the abundance-orbit 
distribution of metal-poor stars in the inner Galaxy, 
Rec <5 kpc. Here, sensible definitions of metal-poor 
could be a) more metal-poor than the oldest disk stars, 
i.e. [M/H] < —1, or b) more metal poor than large pub- 
lished samples of stars with abundances and orbits in 
the inner Galaxy, i.e. [M/H] < —1.5. 

Specifically, we carry out a comprehensive search for 
metal-poor stars towards the Galactic center, drawing 
on the immense wealth of new information that Gaia 
DR3 now affords. Our work relies mainly on newly pub- 
lished low-resolution BP/RP spectra (hereafter collec- 
tively referred to as ‘XP’, Carrasco et al. 2021; De An- 
geli et al. 2022) and photometry synthesized from these 
spectra (Montegriffo, P. et al. 2022a). In principle, var- 
ious metallicity estimates exist for millions of stars in 
DR3. Yet, the XP-based estimates from Andrae et al. 
(2022) suffer from strong biases,” while the RVS-based 
metallicity estimates from Recio-Blanco et al. (2022) are 
currently limited to bright stars with high-quality Gaia 
RVS spectra. 

For our analysis, we select giant stars within 30° of 
the Galactic center, derive data-driven [M/H] estimates 


2 The metallicity biases reported in Andrae et al. (2022) presum- 
ably originate from systematics in the ab initio SED models and 
from our imperfect understanding of the XP instrument used to 
transfer those model SEDs into XP spectra. 


from their XP spectra, and combine them with RVS 
velocities to obtain orbits. This provides a sample of 
1.5 million stars in the inner Galaxy, for a chemody- 
namical study. 

The rest of the paper is organized as follows: in Sec- 
tion 2 we describe the initial selection of the sample, 
and the subsequent data-driven derivation of [M/H] esti- 
mates, using their XP spectra and AllWISE photometry 
(Cutri et al. 2021). In Section 3 we present the spatial, 
[M/H] and orbit distribution of these stars, along with 
the abundance patterns for a small subset. In Section 4 
we summarize these results, put them in the context 
of Galactic archaeology, and sketch a few avenues for 
follow-up work. 


2. DATA SETS AND DERIVED ABUNDANCE AND 
ORBIT QUANTITIES 


We aim to devise a well-defined and large sample of 
stars in the “inner Galaxy” that consists of stellar trac- 
ers that cover the full age and [M/H] range with compa- 
rable selection effects. For that we need tracers that are 
luminous enough to reach the distance of the Galactic 
center and beyond, even in the presence of some dust 
extinction. For chemo-dynamical mapping, we need a 
robust and precise estimate of [M/H] and an estimate 
of the 6D phase space coordinates (xz, v) to estimate the 
orbits. Ideally, we would also like to have precise age 
estimates and individual abundances [X/H], or at least 
an estimate of the a-enhancement, [a/IM/]. Using age es- 
timates and individual abundances for these giant stars 
is beyond the scope of this work, as is a full accounting 
of the sample selection function. 

To build such a sample, we first identify all likely 
red giant stars (RGB and RC) in the direction of the 
Galactic center that: a) are bright enough that we can 
expect a robust and precise estimate of [M/H] from 
Gaia XP data (De Angeli et al. 2022; Montegriffo, P. 
et al. 2022b); and b) have radial velocities from Gaia 
RVS (Katz et al. 2022) to estimate the stars’ orbits in 
conjunction with their parallaxes, w, and proper mo- 
tions, ji. Then we re-derive [M/H] estimates for these 
stars from their XP spectra with a data-driven approach 
that draws on the SDSS’s APOGEE survey (DR17, Ab- 
durro’uf et al. 2022a), to address the documented short- 
comings of the GSP-Phot abundance estimates in Gaia 
DR3 (Andrae et al. 2022). Finally, we calculate the 
orbits of the sample members from their (x, v), which 
requires a model for the Galactic potential, ®(x) (here 
Price-Whelan 2017). 


2.1. Initial Gaia Query 


The initial inner Galaxy sample was devised via the 
following ADQL query: 
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select * from gaiadr3.gaia_source 

where 

parallax < 100.*power(10.,0.2*(0.9 - 
(phot_g_mean_mag - 1.5*(bp_rp-1.)))) 

and parallax < 1. 

and abs(b)<30 and (1<30 or 1 > 330.) 

and bp_rp between 1.0 and 3.5 

and phot_bp_mean_mag < 15.5 

and has_xp_continuous=’ true’ 


The first condition on the parallax selects stars whose 
absolute magnitude is more luminous than Mg = 0.9, 
aimed at selecting giant stars at luminosities that just 
include the red clump (RC); the -1.5* (bp_rp-1.) term 
reflects an approximately dereddened magnitude. The 
query is articulated in such a way that zero or nega- 
tive parallaxes are accommodated. The second paral- 
lax condition only excludes the immediate “foreground” 
at D < 1 kpc. The conditions on / and b select a 
rectangular region of 30° around the Galactic center. 
The lower limit on bp_rp eliminates most luminous hot 
stars for which [M/H]-estimates are not feasible (but 
not all in the presence of reddening). The condition on 
phot_bp_mean_mag is designed to select objects bright 
enough in the blue (BP) part of the XP spectra for ro- 
bust [M/H]-estimates; it still allows to select RC stars 
at the distance of the Galactic center in the presence 
of moderate reddening. The final line in the selection 
ensures that continuous XP spectra are available. This 
query yields 2.1 million objects. 


2.2. Deriving Robust [M/H] Estimates 


While about 75% of this sample has RVS velocities in 
Gaia DR3, only a small fraction of them have metallici- 
ties or abundances derived from the RVS spectra (Recio- 
Blanco et al. 2022). This leaves the route of estimating 
[M/H] from the low-resolution XP spectra of these ob- 
jects. Consequently, our objective is to train a machine- 
learning algorithm that can precisely, accurately, and 
robustly predict [M/H] from Gaia DR3 data. 

Estimating [M/H] based on machine learning involves 
four aspects. 1) It involves the choice of the training 
sample, for which we adopt SDSS DR17 (APOGEE, Ab- 
durro’uf et al. 2022b) since its abundances are well val- 
idated, it covers the inner Galaxy and contains mainly 
giants, many with high extinctions. 2) It involves the 
choice of the data features on which to train the pre- 
diction; these can be derived from the XP spectra but 
may also entail, as we discuss below, external data avail- 
able across the sample, such as near-infrared photometry 
from ALLWISE (Cutri et al. 2021). 3) It involves the 
choice of the machine learning algorithm, for which we 
adopt the extreme gradient boosting algorithm (Chen 
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Figure 1. Top: Validation of our [M/H] estimates from 
Gaia XP data as function of actual [M/H] from APOGEE. 
Bottom: (In-)sensitivity of [M/H] estimates to extinction, 
Ax. The Y-axis shows the difference between [M/H] esti- 
mates based on XP spectra and APOGEE as a function of 
APOGEE’s Ax. There is no evidence for any systematic 
trend of A[M/H] with Ax to an extinction level that corre- 
sponds to Ay 3. 


& Guestrin 2016, hereafter XGBoost). XGBoost is one 
of the best algorithms currently available; it is straight- 
forward and computationally inexpensive to train and 
it can outperform other methods such as deep learning 
(e.g. Grinsztajn et al. 2022). 4) and finally, we need 
to validate our [M/H] prediction, which then speaks to 
the quality of the data, the training sample, and the 
algorithm. 

The choice of the correct data features to train and 
predict [M/H]from the XP spectra is not trivial. One ob- 
vious option is to directly use the 2 x 55 modified Gauss- 
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Hermite coefficients that describe XP spectra. But our 
numerical experimentation implies that other or addi- 
tional data features may yield more precise and robust 
results with a machine learning approach®, for several 
pragmatic reasons. Much of the information about nar- 
row metal line features is contained in high-order co- 
efficients, which are often noisy. We have prior infor- 
mation on which parts of the spectrum are highly diag- 
nostic in [M/H] estimates and which ones are not; we 
need not ask an algorithm to learn this. Any estimate 
of [M/H] requires an implicit estimate of the effective 
temperature Teg, which is covariant with [M/H]; and 
Teg is also highly covariant with reddening. Therefore, 
providing additional information that helps break that 
degeneracy is precious. The near-infrared photometry 
from ALLWISE (Cutri et al. 2021) exists across the sky 
and proves powerful in this respect. 

In practice, we follow the approach of Montegriffo, P. 
et al. (2022a) to calculate synthetic photometry from the 
XP spectra in a wide range of mostly narrowband filters, 
in particular those filters with established utility in iden- 
tifying metal-poor stars: Strémgren filters (Str6mgren 
1966), JPAS+ filters (Marin-Franch et al. 2012) and 
Pristine H& K filters (Starkenburg et al. 2017). On this 
basis, we use XGBoost to train the prediction of [M/H] 
using the entire SDSS DR17 APOGEE of giants with 
Gpp < 15.5 (~ 230,000 stars). Details are given in the 
Appendix. 

Training on a particular data set and using a model 
that is discriminative, rather than generative, raises sev- 
eral issues. First, we (obviously) tie our results to the 
metallicity scale of SDSS DR17. Second, estimates of 
[M/H] may be drawn into the support of the train- 
ing set. In particular, XGBoost will not extrapolate 
[M/H] significantly beyond the APOGEE DR17 range, 
as it is a tree-based method that segments the feature 
space and assigns a mean label to each such segment. 
Finally, generative data-driven spectral models, such as 
The Cannon (Ness et al. 2015), will presumably degrade 
more gracefully towards low signal to noise. Here we 
have addressed the last point implicitly by restricting 
the sample to Gpp < 15.5. 

This results in [M/H] estimates for about 2 million 
giants toward the inner Galaxy, 1.58 million of which 
have RVS velocities, and 1.25 million have both RVS 
velocities and w/o, > 5, which is our minimal condi- 
tion for an orbit estimate. For subsequent analysis, we 
have eliminated 2% of the sample as reddened hot stars, 


3 If stellar features are to be estimated from XP spectra in a full 


forward-modeling approach, the situation may be different. 


which can be readily recognized by their position in the 
m,—bp_rp color plane, where m, = v — 2b+ y is the 
metallicity-sensitive Strémgren filter combination (see 
Stromgren 1966, for details); specifically, we eliminate 
sources with m; + 0.16 (bp_rp— 1) < 0 and bp_rp < 2.7. 


2.3. [M/H] Validation 


We can explore and validate the precision and accu- 
racy of these [M/H] estimates through various compar- 
isons with analogous estimates from (higher-resolution) 
spectroscopic surveys. For this validation, we consider 
SDSS APOGEE (Abdurro’uf et al. 2022a), LAMOST 
(here Xiang et al. 2019), GALAH (DR3 Buder et al. 
2021), and Gaia’s GSP-Spec (Recio-Blanco et al. 2022). 

We start by considering the ~ 17,000 stars in com- 
mon between APOGEE and our sample, shown in Fig- 
ure | (top panel). This panel of Figure 1 illustrates that 
for this sample, XGBoost provides a remarkably precise 
(S 0.1 in the median), accurate and robust [M/H] es- 
timate across —2.5 < [M/H] < 0.5. It should be noted 
that XGBoost was trained on the full all-sky APOGEE 
sample. So this is a validation on a small subset of the 
full training set, not a pure cross-validation. The accu- 
racy of the [M/H] prediction (ie. lack of systematic 
offset) is therefore by construction. It is particularly re- 
markable that basically all stars predicted to be metal- 
poor from [M/H],p are indeed metal-poor according to 
[M/H]lapocgrsr: a selection of, e.g., metal-poor objects 
by [M/H]yp appears to be nearly pure. The bottom 
panel of Figure 1 shows that these [M/H] estimates re- 
main unbiased, at a level of a few percent, even in the 
presence of substantial dust extinction Ax = 0.3, corre- 
sponding to about 3 magnitudes of Ay. 

Appendix A.2 details a true cross-validation of our 
[M/H],p estimates against three external data sets: 
LAMOST, GALAH and GSP-Spec. Broadly speaking, 
the cross-validation affirms that the [M/H]xp estimates 
are precise and robust (enabling pure [M/H]-based sam- 
ple selection). Our estimates have a root mean square 
difference ~ 0.1 dex compared to these datasets, with 
minimal bias. Most importantly, the performance re- 
mains unbiased toward the metal-poor end, affirming 
our confidence in selecting metal-poor stars based on 
XP metallicities. We refer to Appendix A.2 for more 
detailed validation of the [M/H] estimates. 


2.4. Orbits of the Sample Members 


For sample members with suitable 6D phase-space in- 
formation, we calculate their orbits. The sky positions 
and proper motions are available for all of them, radial 
velocities, and good parallax measurements only for a 
good fraction of the sample. In practice, we take all 
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stars that have RVS velocities with dupys < 5 km s7! 


and w/oq > 5, about 1.25 million objects across all 

To compute orbits and compute orbital properties, we 
use a four-component Milky Way mass model consisting 
of spherical Hernquist nucleus and bulge components, 
an (approximate) exponential disk component (Smith 
et al. 2015), and a spherical NFW halo component. We 
adopt a radial scale length hy = 2.6 kpc and scale 
height h, = 300 pc (Bland-Hawthorn & Gerhard 2016), 
and fit for the masses and scale radii of the nucleus, 
bulge, and halo components using the same compilation 
of Milky Way enclosed mass measurements as used to 
define the MilkyWayPotential in gala (Price-Whelan 
2017; Price-Whelan et al. 2020), with the additional con- 
straint of having a circular velocity at the solar position 
Uc(Ro) = 229 km s7! (Eilers et al. 2019). We com- 
pute actions using the “Stackel Fudge” (Binney 2012; 
Sanders 2012) as implemented in galpy (Bovy 2015). To 
compute the orbital eccentricity, pericenter, and apoc- 
enter values, we numerically integrate the orbits with a 
timestep of 0.5 Myr across four times the radial orbit 
period (estimated by the computed orbital frequencies 
from the action solver). 

The [M/H] estimates for 1.5 million stars towards the 
Galactic center, along with orbits for the RVS subsam- 
ple of 1.25 million stars, are available as supplementary 
material to this article and are hosted online’ (Rix et al. 
2022). 


3. RESULTS 


We now present the properties of the low-[M/H] stel- 
lar population contained within this sample. We start 
with the spatial distribution, then move to the [M/H]- 
distribution, the orbit distribution, and finally to a first 
exploration of the abundance patterns of the metal-poor 
stars, especially the level of a-enhancement. This is all 
done with an eye toward what we can learn about a 
central metal-poor in situ halo or bulge population. 


3.1. Spatial Distribution of Metal-Poor Stars 


Figure 2 compares the on-sky density distribution of 
metal-rich stars with [M/H] > —0.4 (left panel) to 
metal-poor ones with [M/H] < —1.5 (middle panel). 
The metal-poor stars exhibit a striking central concen- 
tration relative to the metal-rich sample, suggesting a 
spheroidal population in the inner Galaxy. The pro- 
jected sky density of metal-poor sample members in Fig- 
ure 2 (middle panel) is dramatically altered by the high 
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dust extinction towards the GC (right panel). The fine- 
scale structure in the projected [M/H] < —1.5 stellar 
distribution can mostly be explained by dust extinction, 
as the right-hand panel of Figure 2 shows. This is ex- 
pected as we require Gpp < 15.5 for sample selection. 

We do not have precise parallaxes for the entire sam- 
ple, which complicates the analysis of the spatial struc- 
ture. For the overall sample, the median parallax pre- 
cision is w/o, ~ 12, with 10% of the sample hav- 
ing w/ow S 4, and 1% of the sample with negative 
parallaxes. For the metal-poor sample members with 
[M/H] < —1.5, which turn out to be more distant 
than the sample mean, the median precision is only 
w/Oq~ 5. Therefore, the most robust way to present 
the distance distribution of objects in a directional cone 
may be to consider their parallaxes. Given that the ob- 
jects of particular interest are at ~ 8 kpc distance, it 
matters that we apply a Gaia parallax zero-point cor- 
rection, for which we adopt 0.02 mas (Lindegren et al. 
2021). 

Figure 3 shows n.(a@ | [M/H]) in five [M/H}-bins, with 
the minimal distance (maximal parallax) of 1 kpc at 
the right edge, and the parallax expected for the Galac- 
tic Center indicated by the dotted line (1/8.2kpc, or 
0.12 mas; Gravity Collaboration et al. 2021). This 
figure shows that the distance distribution of giants 
towards the Galactic center depends dramatically on 
[M/H]. Metal-rich stars ([M/H] > —0.4) are distributed 
throughout the disk with heliocentric distances Do = 
1—5 kpc (or 0.2-1.0 mas). In contrast, metal-poor pop- 
ulations become increasingly more concentrated toward 
the Galactic Center (GC). Indeed, the most metal-poor 
stars are predominantly clustered around the GC’s par- 
allax, revealing that there is a centrally concentrated 
very metal-poor population in the Milky Way. 

Figures 2 and 3 qualitatively confirm the picture of a 
metal-poor population that is very concentrated toward 
the Galactic center. The fact that there are almost no 
[M/H] < —1.5 stars at very low latitudes shows that 
almost all lie behind the strong dust extinction 2-3 kpc 
inward of the Sun. Rigorous modelling of these stars’ 
spatial distribution of these stars density distribution 
(e.g. Rix et al. 2021) is beyond the scope of this paper, 
with the severe dust extinction and crowding effects pre- 
senting formidable challenges. However, we can model 
the parallax distribution (Fig. 3) in the central regions 
less affected by dust to quantify the radial extent of this 
metal-poor population. 

We assume that the 3D distribution is a spheri- 
cal Gaussian as a function of Galactocentric radius R 
PGauss(R, ORgc), With a width of or,,. As long as we 
can see to the far side of the Galactic center (see Fig. 3), 
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Figure 2. On-sky logarithmic density distribution of stars in our sample for metal-rich (left panel) and metal-poor (middle 
panel) selections. The distribution of metal-poor stars appears very centrally concentrated, but the morphology is dramatically 
modulated by the foreground dust extinction, as illustrated in the right panel. The right panel demonstrates that the dearth of 


stars is highly correlated with the distribution of dust. 
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Figure 3. Parallax distribution of the different sub-samples, 
selected by [M/H], taking w as a proxy for distance. It is 
apparent that the metal-poor samples are centered around 
the expected parallax of the Galactic center (dotted black 
line; Gravity Collaboration et al. 2021), ie. they are very 
concentrated near the Galactic Center (GC), increasingly so 
towards lower [M/H]. 


we should expect the parallax (@) distribution in a cone 
of dQ in the direction (J, b) to be 


n(w) = dNdww* poauss(R(w,1,b),oRec)- (1) 


Phrasing this model in terms of the @ allows the model 
to deal with vanishing (or even negative) observed par- 
allaxes. We find that the maximum likelihood of the 
observed parallaxes (and their uncertainties) for the 


[M/H] < —1.5 stars implies an extent of org, ~ 2.7 kpc. 
Interestingly, this radius corresponds to an angle of 18° 
at the distance of the Galactic center, which seems very 
plausible in light of the projected density distribution 
(Figure 2; middle panel). 


3.2. The [M/H]-Distribution in the Inner Galaxy 


We now turn our attention to the [M/H]-distribution, 
nx([M/H]), in the inner Galaxy, which is illustrated in 
Figure 4. The first remarkable feature of the distri- 
bution is the sheer number of objects, most apparent 
from the cumulative distribution in Figure 4: there are 
> 4,000 stars with an estimated [M/H] < —2, ~ 18,000 
with [M/H] < —1.5, and almost 100,000 stars with 
[M/H] < —1, ie. below the abundance floor of the old, 
a-enhanced disk. This is an order of magnitude more 
objects than previously published metal-poor samples 
([M/H] < —1) of the inner Galaxy (e.g., Arentsen et al. 
2020a,b). 

A rough estimate shows that this corresponds to = 
5 x 10’Mo in stars at [M/H] < —1.5 in the inner few 
kpc of the Galaxy. Taking the nearest globular cluster 
M4 as a template for an old, [M/H] < —1 population, 
we find that there are about 35 giants with Mg < 0.5 
in the Gaia DR3 catalog beyond M4’s half mass ra- 
dius of 6’ (Richer et al. 2004, where the Gaia catalog 
should be approximately complete). Given M4’s to- 
tal mass of 7 x 104Mo (Marks & Kroupa 2010), this 
implies one giant (at Mc < 0.5) per 800Mo of to- 
tal stellar population mass. At face value, this implies 
Mtot((M/H] < -1.5) ~ 2x 10’Mo. But this does not 
account for the incompleteness of XP spectra in very 
crowded fields, not for the far-reaching effects of dust 
extinction. A conservative estimate of these effects is a 
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factor of 2.5, leading us to a lower limit on the mass of 
Mtot([M/H] < —1.5) = 5 x 107Mo. 

We have already presented two arguments that the 
low-[M/H] portion of n.([M/H]) is not severely contam- 
inated by spurious [M/H] estimates from some of the 
far more numerous stars of higher (true) metallicity: 
the validation with APOGEE and the parallax, or dis- 
tance, distribution, which is most peaked for the most 
metal-poor objects. The [M/H] distribution itself pro- 
vides additional evidence. It rises very steeply from 
[M/H] = —2.5 to —2.2. The distribution then follows 
a power law of slope d(logn.)/d[M/H] = 1 over a wide 
range of metallicities to [M/H] ~ —0.9, where it steep- 
ens towards a peak at [M/H] = —0.6, beyond which 
nx([M/H]) starts dropping towards the highest metal- 
licities present, [M/H] ~* 0.5. 

The bottom panel of Figure 4 shows the metallicity 
distribution with stars of an estimated Galactocentric 
distance Rac < 4 kpc in green, and the stars of the 
intervening disk in gray. As giants, in particular red 
clump stars, are of comparable luminosity and color over 
a wide range of metallicities, the differential selection 
effects across different [M/H] should be modest. Hence, 
the [M/H]-distribution of stars Rec < 4 kpc in Figure 4 
should be relatively unbiased. The [M/H] distribution 
of course varies with position, in the Galaxy (both |z| 
and Rac), as Figure 4 illustrates. 

The decrease in n.([M/H]) for stars at the metal-rich 
end ({M/H] > 0) with Rac < 4 kpc (green) may seem 
surprising, because we know from APOGEE (e.g. Eilers 
et al. 2022) that the inner Galaxy is teeming with high- 
metallicity stars. This decrease in n,.([M/H)]) is straight- 
forwardly explained by Figure 2, which shows that in the 
inner Galaxy, we are almost completely lacking sample 
members with |Z| < 700 pc; any thinner population 
will inevitably be absent from our sample due to dust 
extinction. The APOGEE sample has a different MDF 
in part because its 2MASS-selected stars better sample 
the low-|Z|, high-[M/H] population. 

The bulk of the metal-poor population, below the min- 
imal [M/H] of the old disk at [M/H] ~ —1, follows 
d(log n,.)/d[{M/H] = 1: the cumulative number of stars 
N,.(< Z) below a (linear) metallicity Z grows linearly 
with Z. This seems to be a natural slope for n,.([M/H]) 
in the early phases of a self-enriching, in situ system, as 
shown by the following argument, which builds on the 
models of Weinberg et al. (2017). 

Under fairly general conditions, the evolution of the 
ISM metallicity in a simple, fully mixed system can be 
described by 


dMz 
dt 


= UZ M, —(l+n-r)M,.Z, (2) 


10° 


— 
So 
or 


eK 
j=) 
is 


# of Stars 
2, 


-1 
[M/H]xp 

Tt 1 

Ree < 4kpc 

Rec > 4kpc 


# of Stars 


— 
j=) 
NO 

T T Ported 

i i rill 


— 
j=) 
an 


-2 -1 


[M/H]xp 


Figure 4. Top: Cumulative distribution of metallicities 
[M/H] for the entire sample of ~ 1.7 million objects, es- 
timated from XP spectra and broad-band photometry via 
XGBoost prediction trained on the APOGEE DRI17: there 
are ~ 5,000 stars with [M/H] < —2, ~ 20,000 stars with 
[M/H] < —1.5, and nearly 100,000 stars with [M/H] < —1. 
Bottom: Distribution of [M/H] for our sample, split into 
sources with presumed Galactocentric radii Rac < 4 kpc 
(in green), which are the focus of our analysis, and those at 
larger Rac. The (green) distribution shows a steep rise to 
[M/H] ~ —2, then a power-law of slope d(logn.)/d{[M/H] ~ 
1 to [M/H] ~ —1, and finally a yet again steeper rise to 
[M/H] ~ —0.6. A possible interpretation of these slopes in 
the context of a simple chemical evolution model is given in 
Section 3.2. 
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where Mz is the mass in metals, yz is the IMF-averaged 
yield, M, is the star formation rate (SFR), Z = Mz/M, 
is the metal mass fraction of the star-forming gas, 7 = 
Mout / M, is the mass loading factor of a gaseous outflow, 
and r 0.4 represents the recycling of metals that are 
formed into stars but quickly returned when the stars 
die. Equation (2) treats enrichment as instantaneous, 
which should be a good approximation at low metallic- 
ities where core collapse supernovae and massive star 
winds are the dominant sources. It also assumes that 
outflows have the same metallicity as the ambient ISM; 
if outflows instead consist of a fraction f of the super- 
nova ejecta plus entrained ISM, then the effect is to 
reduce yz by a factor (1 — f) so that it reflects only 
metals retained by the ISM. Importantly for the case 
at hand, if Z is well below the equilibrium metallicity 
Zeq = yz/(1+n-— 1) then the entire second term of 
equation (2) can be neglected. 

In this low metallicity limit, equation (2) is simply 
Mz = yzM,, and for a metallicity-independent yield, 
its time integral implies Mz = yzM, and thus Z = 
yzM,,/M,. The log-slope of the MDF is 


dlog M, Z My YZ 


= = = 3 
dlog Z M, Z Zr. (3) 


dlog nx 

d{M/H] 

where the last equality introduces the star formation 
efficiency (SFE) timescale 7, = M,/M.. Using 
d M, M.M, 


Z= ay (Mz /Mz) vay ~ v2 ME (4) 


gives the end result 


dlog nx 7 Zt ; x, Ie (7) TM, 
d(M/H] YZ M, M, 

| (5) 
The combination (7.M,/M,) corresponds to the frac- 
tional change of the gas reservoir mass over one SFE 
timescale. 

Equation (5) shows than an MDF log-slope © 1 is a 
generic result in the low metallicity regime (Z < Zeq), 
arising whenever the gas mass is constant (M, = 0) 
or more generally when the star-to-gas ratio M,/M, is 
small enough that the second term in Equation (4) can 
be neglected. Figure 4 exhibits this generic slope over 
the range —2 < [M/H] < —1. The steeper slope at 
[M/H] < —2 could plausibly arise because the gas reser- 
voir is small but rapidly growing, with 7,M,/M, > 1. 
The steeper slope at —1 < [M/H] < —0.6 may re- 
flect rapid gas accretion (high M,/M,) coinciding with 
the onset of the old a-enhanced disk (e.g. Belokurov & 
Kravtsov 2022; Xiang & Rix 2022; Conroy et al. 2022). 
At still higher metallicities our neglect of sink terms in 
equation (2) becomes a poor approximation. 


3.3. The Orbit Distribution of low-[M/H] Stars in the 
Inner Galaxy 


We now turn to the orbit distribution of the sample 
at hand (see Section 2.4), which can tell us which of the 
metal-poor stars selected in the inner Galaxy remain 
confined to the central regions of the Milky Way. The 
orbit distribution can tell us which stars are just passers- 
through near their pericenter, but on orbits that take 
them into the outer halo. In particular, we might ex- 
pect to see the pericenter members of accreted halo com- 
ponents like GSE, with stars on highly eccentric orbits 
that take them to Rapo > 10 kpc. We can also delineate 
how early — or how metal-poor, if we use abundances 
to estimate relative stellar ages — the population devel- 
ops net rotation and how rapidly the kinematics change 
with [M/H] towards near-circular motion in the plane 
of the Galaxy (see Arentsen et al. 2020a; Belokurov & 
Kravtsov 2022; Conroy et al. 2022). 

We start by considering the orbit distribution, 
nx(e, Rapo), in terms of the orbital eccentricity, e and the 
apocenter Rapo. We initially focus on the [M/H]-range 
—1.9 < [M/H] < —1.2 since this encompasses the bulk 
of known GSE stars. This distribution of 28,000 stars is 
shown in Figure 5. The density in the (e, Rapo)-plane is 
affected and limited by two experimental aspects. First, 
our initial Gaia query selects almost exclusively stars 
with Rac < 8 kpc, and hence all stars whose pericen- 
ter distance exceeds 8 kpc are excluded; this boundary 
is indicated by the thick gray line. Except at e < 0.2 
the n.(e, Rapo) distribution falls off toward large Rapo 
well before this limit. There may be two reasons for 
the seeming dearth of stars with Rap. < 2 kpc near 
the center. First, severe dust extinction simply obscures 
much of the central kiloparsec (Fig. 2). Second, any 
uncertainties in the 6D phases-space coordinates, espe- 
cially uncertainties in the distance, are far more likely 
to increase the inferred Rapo than decrease it, as it is a 
positive definite quantity. These aspects deserve careful 
modeling, which is beyond the scope of this paper. 

The distribution in Figure 5 can be characterized by 
two regimes: most stars form a nearly flat distribution 
in eccentricity from e = 0.1 to e = 0.8, but with a 
narrow range in Rapo, with (Rapo) * 4 kpc. This im- 
plies that most stars form an approximately isotropic 
distribution that stays confined to the inner Galaxy and 
hence will not appear in surveys of halo stars that fo- 
cussed on larger Rac and off the Galactic plane. The 
second regime is that of e > 0.75 and Rapo 2 10 kpc. 
The orbits of these stars fully match the expectations 
of the pericenter members of the GSE stars, the ac- 
creted component that dominates the halo population 
between 10 and 30 kpc (Helmi et al. 2018; Belokuroyv 
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Figure 5. Top: Rapo vs. eccentricity space colored by the 
median XP metallicity. We overlay contours of linear density 
using a Gaussian kernel density estimate. Bottom: Distri- 
bution of the —1.9 < [M/H] < —1.2 sub-sample in Rapo vs. 
eccentricity space. This figure shows that most of the low- 
metallicity sample members have orbits that remain confined 
to the inner Galaxy (Rac < 5 kpc), with a broad eccentric- 
ity distribution. At high eccentricities, there are many stars 
with large apocenters, potentially representing a population 
of accreted stars on highly radial orbits. 


et al. 2018a; Naidu et al. 2021). Indeed, recent work 
by Belokurov et al. (2022) has shown quite clearly that 
there are likely GSE members with pericenters within 
3 kpc. There is tantalizing but inconclusive evidence 
from Figure 5 that the distribution of GSE stars might 
extend to small Rapo and lower e: Figure 5 indicates 
an excess of stars near (e, Rapo) © (0.7,4 kpc). Here it 
would become important to disentangle GSE from other 


proposed accreted structures in the inner Galaxy (e.g., 
Kruijssen et al. 2019; Horta et al. 2021). 

We now ask to what extent different metal-poor stars 
in the inner Galaxy have net angular momentum. In 
general, more metal-rich stars are generally expected 
and observed to have a more coherent sense of angu- 
lar momentum. We explore this in the two panels Fig- 
ure 6, which show how close, or far, stars at a given 
[M/H] and Rapo are from being an ensemble of stars on 
circular, in-plane orbits. This is quantified via the ra- 
tio of the angular momentum, which is the azimuthal 


action Jy, to the total action Sto, = \/J5 + JR + J?. 


The top panel shows that the increase of Jg/Jtot with 
[M/H] depends somewhat on the orbit size: At a given 
[M/H]orbits with smaller R,,. are farther away from the 
coplanar quasi-circular orbits than at Rapo ~ 5 kpc. 

The bottom panel of Figure 6 shows J¢/Jtot({[M/H)]), 
marginalized over 3 < Rapo/kpe < 7: starting al- 
ready at [M/H] < —2 there is some, albeit small, net 
J¢/Jtor({M/H]), which increases gradually to ~ 0.3 
at [M/H] = —1.3. Among more metal-rich popula- 
tions, J/Jtot([M/H]) steeply rise with rising [M/H] to 
[M/H] ~ —0.7, after which it levels off in the regime 
dominated by cold rotation. This analysis implies that 
some modest net rotation is present in the popula- 
tion well below [M/H] < —1.5, affirming and extend- 
ing the findings of Arentsen et al. (2020a), Belokurov & 
Kravtsov (2022) and Conroy et al. (2022). 


3.4. Abundance Patterns in the Inner Galaxy 


We now enlist element abundances for our sample 
stars that overlap with APOGEE to explore which parts 
of the sample appear accreted from a chemical enrich- 
ment perspective (as opposed to proto-Galactic, or in 
situ). Such abundance-based arguments for the origin of 
stars are based on the idea that enrichment works dif- 
ferently in (sub-)halos of different mass (e.g., Hawkins 
et al. 2015; Das et al. 2020; Horta et al. 2021; Belokurov 
& Kravtsov 2022): the production of a-elements and Al 
aided by high star-formation intensities in the deeper po- 
tential wells of the proto-Galaxy, compared to the shal- 
lower potential of lower-mass satellite galaxies. How- 
ever, these abundance diagnostics of the stars’ forma- 
tion environment may be less discriminating at lower 
metallicities (e.g., Conroy et al. 2022). 

We remove stars belonging to globular clusters and ap- 
ply a few basic quality cuts to the APOGEE abundance 
data based on their uncertainties and flags. We plot 
the resulting chemical distributions, [X/Fe] vs. [Fe/H] 
with X= {a, Al, Mn}, as a function of eccentricity in 
Figure 7, , by median apocenter radius. The presumed 
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Figure 6. Dependence of the orbits’ level of rotation sup- 
port, as a function of [M/H] and Rapo. The rotation support 
is defined as (J¢/Jtot), where 1 is an ensemble of circular, 
in-plane orbits, covering both an increase in net rotation 
and the approach towards circular, co-planar orbits. The 
top panel shows (J%/Jtot) as a function of [M/H] and Rapo. 
The bottom panel shows J¢/Jtot({[M/H]), marginalized over 
3 kpe < Rapo < 7 kpe: starting already at [M/H] < —2 
there is some net J¢/Jtot([M/H]), which increases to ~ 0.3 
at [M/H] = —1.3, then steeply towards [M/H] ~ —0.7, after 
which it levels off in the cold rotation dominated regime. 


boundaries between in situ. and accreted are indicated 
by the yellow dashed lines. 

Figure 7 shows clearly that at modest eccentricities 
(left and middle columns), most stars belong to the a- 
enhanced sequence, implying a proto-Galactic origin. At 
high eccentricities (right columns), a prominent low-a 
sequence with large apocenters appears. These stars 
on eccentric orbits with large apocenter appear to be 
remnants of the accreted GSE merger based on their 
chemistry and kinematics (e.g., Mackereth et al. 2019; 
Naidu et al. 2020; Bonaca et al. 2020; Hasselquist et al. 
2021; Horta et al. 2022). Note that among the stars of 
highly eccentric orbits there is a portion whose abun- 
dances point towards proto-Galactic origin (right col- 
umn). Remarkably, they almost exclusively have apoc- 
enters < 5 kpc, very much like their cousins on less 
eccentric orbits, but very unlike the chemically iden- 
tified GSE debris. The population of eccentric stars 
attributable to the proto-Galaxy has been seen before. 
It has been dubbed ‘Aurora’ by Belokurov & Kravtsov 
(2022) and has been explored at even lower metallicities 
by Conroy et al. (2022). Our sample probes to consider- 
ably lower Rac than either of these works, and reveals 
the bulk and full extent of the metal-poor proto-Galaxy 
at the heart of the Milky Way. 

It is apparent from Figures 5 and 7 that both orbits 
and abundances can help to differentiate between pre- 
sumed proto-Galactic stars and stars accreted later from 
a lower-mass satellite. A differentiation by chemistry 
alone appears increasingly difficult at low [M/H] (Con- 
roy et al. 2022). But the combination of both aspects, 
summarized in Figure 7, can make such a differentiation 
more convincing. 

Horta et al. (2021) selected “accreted” stars in these 
APOGEE chemical spaces as evidence of a past ‘Hera- 
cles’ merger buried deep within the MW, distinguished 
from GSE on the basis of a bimodal distribution in to- 
tal orbital energy. This claimed galaxy bears resem- 
blance to the previously proposed Kraken/Koala merg- 
ers (Kruijssen et al. 2019, 2020; Forbes 2020), and these 
names may well refer to the same event. However, Lane 
et al. (2022) pointed out that the apparent bimodality 
in orbital energy used by Horta et al. (2021) to select 
Heracles/Kraken/Koala might be an artifact of the spa- 
tial selection function of APOGEE. Furthermore, Her- 
acles/Kraken/Koala and the proto-Galactic Aurora are 
virtually indistinguishable in all chemical spaces (Horta 
et al. 2022; Naidu et al. 2022; Myeong et al. 2022). 

Taken together, these lines of evidence suggest that 
the stars attributed to Heracles, Kraken, Koala and Au- 
rora, along with the numerous metal-poor stars iden- 
tified here are different orbit and metallicity regimes 
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Figure 7. Abundance diagnostics for stars in the inner Galaxy, and their relation to these stars’ orbital apocenters. We show 
[Fe/H], [a/Fe] (top row) and [Mg/Mn] and [Al/Fe] (bottom row) from SDSS/APOGEE (DR17) for stars with @ < 0.25 mas. 
The three columns of panels show slices of increasing orbital eccentricity from left to right, with stars colored by their median 
Rapo. The color map transitions from blue to red for apocenters at the Solar Galactocentric radius; linear density contours of the 
underlying sample are overlaid in black. The gold dashed lines indicate boundaries that have been used for selecting ‘accreted’ 
stars on the basis of their abundance patterns (e.g., Hawkins et al. 2015; Das et al. 2020; Horta et al. 2021; Belokurov & 
Kravtsov 2022; Conroy et al. 2022). Note, however, that the [Mg/Mn]-[Al/Fe] diagnostic has been found to be more ambiguous 
at low metallicity (see Belokurov & Kravtsov 2022; Conroy et al. 2022). Overall, the figure shows that there is a high degree of 
correlation between the orbital and the abundance properties. At eccentricities e < 0.8 most stars with [M/H] < —1 have both 
Rapo < Rac(Sun) and abundance patterns attributed to in situ formation. Only at the highest eccentricities, e > 0.8, is there 
a good fraction of stars that have large apocenters and at the same time abundance patterns characteristic of accreted stars. 
Both of these latter properties are expected for members of GSE near pericenter. 
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of one proto-Galactic stellar population, now forming 
the metal-poor, old heart of the Milky Way. Whereas 
past studies have mainly identified the eccentric tail of 
this population that traverses the solar neighborhood 
or off the Galactic plane (Belokurov & Kravtsov 2022; 
Conroy et al. 2022; Myeong et al. 2022), we now show 
more directly that the bulk of this population resides in 
the bulge itself, mostly within a few kpc of the Galac- 
tic center. Indeed, this matches theoretical predictions 
from cosmological zoom-in simulations, which predict 
a spheroidal distribution of proto-Galazy stars predom- 
inantly concentrated around the Galactic center (e.g., 
El-Badry et al. 2018, Wetzel et al. 2022, Figure 13 in 
Belokurov & Kravtsov 2022). 

The picture presented here appears to obviate the 
need for an ex-situ merger like Heracles/Kraken/Koala, 
at least based on the orbits and chemistry of field stars. 
This does not necessarily imply that these mergers did 
not occur; rather, it emphasizes the — conceptual and 
practical — ambiguity between in situ and accreted stars 
during the earliest phases of the Milky Way, which led 
us to use the collective term proto-Galaxy. This view is 
supported by cosmological zoom-in simulations, which 
affirm the difficulty in distinguishing accreted galaxies 
from the MW proto-galaxy at early times (e.g., Renaud 
et al. 2021a,b; Orkney et al. 2022). 


4. SUMMARY, DISCUSSION AND OUTLOOK 


We have presented an exploration of the metal-poor 
stellar population in the inner Galaxy (Rac < 5 kpc), 
drawing on the newly available low-resolution XP spec- 
troscopy from Gaia’s DR3. The primary goal was to 
see whether there is an extensive ancient and metal- 
poor stellar population at the “heart” of the Milky Way. 
We wanted to learn whether the bulk of these stars are 
a tightly bound proto-Galactic population, rather than 
stars accreted from more distant satellites, whose lo- 
cation in the inner Galaxy just reflects the pericenter 
phase of orbits that take them out to the “classical” 
halo at Rec 2 Ro. As simulations imply that the ear- 
liest heart of galaxies often arises from the high-redshift 
coalescence of clumps with comparable mass, we pre- 
fer the term proto-Galactic for the resulting population, 
rather than in situ, a term that singles out one of these 
clumps even if it has only modestly larger mass than 
other nearby clumps. 

The practical foundations of our analysis were esti- 
mates of [M/H] of ~ 2 million giant stars within 30° 
of the Galactic center, at least 1 kpc from us and suf- 
ficiently bright for good S/N in the blue part of the 
XP spectra, Gpp < 15.5. We converted the XP spec- 
tral information into a set of narrow-band fluxes in fil- 


ters that are known to be highly diagnostic of [M/H] in 
cool stars. We combined this synthetic photometry with 
ALLWISE photometry and then used them — along with 
the XP spectral coefficients — for a data-driven [M/H] 
prediction, trained on the APOGEE DR17 data using 
XGBoost. We found these [M/H] estimates to be pre- 
cise to ~0.1 dex. We also found them robust and pure 
enough to identify metal-poor samples of stars in the 
inner Galaxy ([{M/H] < —1), despite the presence of a 
vastly dominant population of more metal-rich stars in 
the inner Galaxy. For most of these stars, radial veloci- 
ties exist from Gaia RVS, allowing the estimates of their 
orbits. 

On this basis, we identified metal-poor samples of 
stars towards the inner Galaxy that are about two orders 
of magnitude larger than published ones: > 4,000 stars 
with [M/H] < —2, ~ 20,000 stars with [M/H] < —1.5, 
and ~ 70,000 stars with [M/H] < —1. The most 
metal-poor tail of the metallicity distribution extends 
to [M/H] » —2.5. Whether this constitutes a genuine 
metallicity floor or reflects that the APOGEE training 
set for XGBoost only extended to [M/H] = —2.5 requires 
follow-up data. 

These samples reveal a large metal-poor population 
({M/H] < —1.5) in the inner few kpc of our Galaxy. 
Three lines of argument imply that this population is 
centrally concentrated: first, the parallax distribution 
implies that most stars are within 5 kpc of the Galac- 
tic center (Figure 3). If we model the spatial density 
distribution of stars with [M/H] < —1.5 as a Gaussian 
centered on the Galactic center, with an extent of oRQo, 
this parallax distribution implies org, ~ 2.7 kpc. Sec- 
ond, the number density projected onto the sky shows a 
centrally concentrated distribution (albeit severely mod- 
ulated by intervening dust, Figure 2), which is qualita- 
tively consistent with org, ~ 2.7 kpc. Third, their 
orbits show that most stars have apocenters of less than 
5 kpc (Figure 5). The latter property implies that most 
of this population could not have been found in many 
of the previous surveys that focused on Rac > 5 kpc. 
At the same time, our analysis showed that a minority 
(but still a large number) of stars with [M/H] < —1.2 
currently in the inner Galaxy are on highly eccentric 
orbits that can take them to > 10 kpc. These are pre- 
sumably members of the accreted halo, found here near 
the pericenter of their orbits. 

Using the members of this sample that have de- 
tailed abundances from SDSS (DR17), in particular 
[a/Fe], [Al/Fe] and [Mn/Fe], we explored which of these 
stars have abundance patterns attributable to a proto- 
Galactic (or, in situ) vs. accreted origin. We found that 
almost all stars that remain tightly bound within the 
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inner Galaxy (Rapo < Rac(Sun)) have abundances that 
point towards a proto-Galactic origin. Stars on highly 
eccentric orbits with Rapo > 10 kpc show abundances 
typical of accreted stars and are presumed pericenter 
members of GSE. 

We quantified how much net rotation the various 
mono-abundance populations show and how much they 
resemble disk-like orbits. We did this by considering 
the mean ratio of the angular momentum to the total 
action, Jy/Jtot([M/H]). We found that only the most 
metal-poor tail [M/H] < —2 of this population shows 
no rotation, Jy/Jtot([M/H] < —2) =» 0. Populations 
with higher [M/H] also have higher Jg/Jto:([M/H]): at 
[M/H] = —1.5, Je/Jtot has reached 0.25, then continu- 
ously increasing to J/Jtot ¥ 0.9 in [M/H] = —0.8. 

All of this information fits a picture in which this 
metal-poor heart of the Milky Way constitutes the most 
ancient proto-Galactic component of our Galaxy. 


e Much of this population is confined to Rac < Ro, 
or tightly bound at the bottom of our Galaxy’s 
potential well, arguing against an origin through 
accretion from a once-distant satellite. 


e The vast majority of stars with [M/H] < —1.5 
are distinctly a-enhanced, arguing for rapid en- 
richment in a deep potential well. 


e The stellar mass we see at [M/H] < —1.5 within 
5 kpe corresponds to ~ 5 x 10’Mo, but manifest 
dust obscuration implies a large correction that 
suggests that the stellar mass of this metal-poor 
heart of the Milky Way is high, M, > 10°Mo. 
Such high masses of strongly a-enhanced stars at 
such low [M/H] are most likely to occur at the 
center of a (eventually) quite massive halo, not at 
a low-mass satellite. 


e We have mapped the degree of (prograde) ordered 
azimuthal motion, J¢/Jtot([M/H]) as function of 
[M/H]. At [M/H] ~ —2 it is approximately zero, 
and increases continuously when considering sets 
of stars with increasing [M/H]. Such a relation is 
expected if star formation occurs — as time goes on 
—in gas that has increasingly settled into a disk, as 
recent simulations of disk galaxy formation show 
(Gurvich et al. 2022). 


e Much of this central component has [M/H] well 
below that of the centrally concentrated, a- 
enhanced, thick disk, whose oldest members (at 
[M/H] ~ —1) appear to be = 12.5 Gyrs old (Xi- 
ang & Rix 2022). If much of this metal-poor heart 
of the Milky Way chemically, and therefore tempo- 
rally, predates the oldest a-enhanced disk, it must 


have formed at 7 2 12.5 Gyrs, corresponding to 
z 25. This stellar population at the heart of the 
Milky Way should be almost exclusively ancient. 


The results presented here are by no means a new dis- 
tinct stellar component of the Milky Way. The stars in 
our sample seem to constitute the tightly bound part — 
and bulk — of a proto-Galactic spheroid. The distribu- 
tion tail of stars on somewhat more extended orbits has 
already been recognized in recent work as in situ halo 
(Belokurov & Kravtsov 2022; Conroy et al. 2022). But 
our results significantly flesh out the existing picture by 
showing that there is indeed a tightly bound in situ “ice- 
berg”, whose tips have been recognized before. The fact 
that there are many metal-poor stars in the inner galaxy 
has also already been recognized by the recent work of 
Arentsen et al. (2020a,b). 

Our analysis here is in many ways preliminary and 
suggests various avenues of follow-up. The spatial distri- 
bution of this population deserves to be modeled quan- 
titatively, including the effects of dust extinction, source 
crowding, parallax uncertainties, and the giant luminos- 
ity function. The abundances, and perhaps more im- 
portantly, the abundance patterns, deserve full spectro- 
scopic follow-up to understand: How pure are the XP 
selected samples in this regime? Is there a floor in the 
[M/H] distribution that might reflect the [M/H] of the 
proto-Milky Way’s circumgalactic medium? If these are 
indeed among the most ancient stars in the Milky Way, 
do they have remarkable abundance ratios? Was the for- 
mation of this component associated with the formation 
of an extensive set of globular clusters? What can the 
degree of central concentration of this component tell us 
about the violence of subsequent merging events, which 
would scatter stars to larger radii? 

On a more technical side, this analysis reflects the 
astounding information content of the Gaia DR3 data, 
particularly the XP spectra. Our data-driven approach 
to estimate [M/H] seems to work well for the current 
analysis, at the price of restricting [M/H] estimates to 
bright objects. Presumably, we are far from exploiting 
the information of the XP spectra, which should be un- 
locked by forward-modeling of the data. 
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APPENDIX 


A. [M/H] ESTIMATES VIA XGBOOST TRAINED 
ON SDSS APOGEE ABUNDANCES 


Here we expand on a few specifics of the [M/H] esti- 
mate that is summarized in Section 3.2. 


A.l. Training of XGBoost model of [M/H] 


We want to derive [M/H]from XP information for ob- 
jects that still have significant flux in the blue (Ggp < 
15.5), yet may be reddened by several magnitudes. To 
help break any degeneracies between extinction and Teg 


we include near-infrared photometry from ALLWISE 
(Cutri et al. 2021). We adopt SDSS DR17 APOGEE 
abundances (Abdurro’uf et al. 2022b) as a training sam- 
ple, since it covers the inner disk and contains many 
giants stars with high extinctions. 

For the machine learning algorithm, we choose the 
extreme gradient boosting algorithm (Chen & Guestrin 
2016, hereafter XGBoost), as it is computationally inex- 
pensive to train and can outperform other algorithms, 
e.g. Deep Learning (e.g. Grinsztajn et al. 2022). To im- 
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Table 1. Test errors of [M/H] from XGBoost from 20-fold 
cross-validation using various numbers of XP coefficients. 
We quote the median absolute error and the root-mean- 
square error. 


input features median AE RMSE 


all 55 coefficients 0.061 0.145 
first 40 coefficients only 0.059 0.138 
first 30 coefficients only 0.060 0.137 
first 20 coefficients only 0.062 0.139 
first 10 coefficients only 0.066 0.144 
XP colours only 0.078 0.154 
all coefficients + XP colours 0.064 0.138 
all coefficients + WISE 0.088 0.186 
XP colours + WISE 0.067 0.136 
all coefficients + XP colours + WISE 0.049 0.107 


prove the [M/H] predictions for our sample (consisting 
of red giant stars by construction), we also restrict the 
training sample from SDSS DR17 sample to only giants, 
via log g < 3.5. 

This leaves us with a suitable choice of data features, 
for which we explored a number of options. Our start- 
ing point was the set of XP coefficients (De Angeli et al. 
2022; Carrasco et al. 2021), which we normalized by 
the Gaia apparent G flux. In Table 1, we quote the 
test errors for using the different numbers of XP coeffi- 
cients. The performance is overall very good but depend 
only weakly on how many coefficients we use. This sug- 
gests that also the high-order coefficients contain useful 
information, and, at least for this application and the 
APOGEE DR17 sample, a truncation of XP coefficients 
is neither required nor helpful in our case. 

Given the results in Montegriffo, P. et al. (2022a), we 
then attempted to only use photometry synthesized from 
XP spectra using GaiaXPy° and combine this with ALL- 
WISE near-infrared photometry. Specifically, we syn- 
thesize photometry from the systems of Pristine, Strom- 
grenStd, Jpas, Jplus, SkyMapper, and HstAcswfc. As 
we can see from Table 1, using only such colors has a 
similar performance as only using the XP coefficients, 
despite including ALLWISE photometry. However, a 
substantial improvement is obtained when combining all 
XP coefficients with these colors. 

In the end, we used GaiaXPy to compute the follow- 
ing synthetic photometry and combined it with all 110 
XP coefficients in the XGBoost training on SDSS DR17 
(APOGEE): 


5 https: //gaia-dpci.github.io/GaiaX Py-website/ 


G-W, 
Gpp — Grp 
Gpp — W2 
W, —- W2 


CaHK pristine _ Wi 

CaHK pristine — UStromgrenStd 

CaHK pristine _ bstromgrenStd 

CaHK pristine — YStromgrenStd 

Jplus, — Jplus; 

Jplusy395, — CaHK pristine 

Jplusys;5 — CaHbK pristine 

Jplusggg; — CaHKpristine 

Grp — Jplus, 

UStromgrenStd — 2: bstromgrenStd + YStromgrenStd 
CaHK pristine — 2° Jplusg419 + Jplusys3q 


We emphasize that the Ggp — W2 color has the longest 
wavelength leverage and therefore is highly affected by 
extinction, whereas the W, — W2 color is comparatively 
insensitive to extinction. 

As is evident from Table 1, this final configuration has 
slightly worse performance than the best configuration. 
However, the median absolute errors are similar, that 
is, most sources get similar results. The only notewor- 
thy difference is in the RMS error, suggesting that the 
final configuration may have a few more outliers. Still, 
the results are better than what we can achieve from 
coefficients alone. 


A.2. Further [M/H] Validation 


Figure 8 illustrates a cross-validation comparison of 
our results to other spectroscopic surveys. Such cross- 
validation is easy to interpret if one can assume that 
the external data set represents “ground truth”. We ap- 
plied a SNR_G > 25 to the LAMOST sample and required 
that there were no quality (concern) flags in any of the 
GALAH spectra. The left panel shows that the [M/H]xp 
estimates are accurate for LAMOST, which itself is tied 
to the APOGEE [M/H] scale, with a slight increase 
in scatter towards low [M/H]. The comparison with 
GALAH (middle panel) affirms the precision and pu- 
rity of the low-[M/H]xp estimates. Note that GALAH 
has a different scaling between [M/H] and [Fe/H] than 
APOGEE and LAMOST, which may explain the offset 
of the estimates for low-[M/H] stars from the 1-to-1 line. 
In the parent sample without stringent GALAH quality 
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Figure 8. Cross-validation of our [M/H]xp estimates against three external spectroscopic data sets, LAMOST (here Xiang 
et al. 2019) , GALAH (Buder et al. 2021) and GSP-Spec (Recio-Blanco et al. 2022). This cross-validation supports the precision 


and robustness of our [M/H]xp estimates. 


cuts, there are many stars for which GALAH implies 
that they are [M/H] < —1, but the [M/H]xp estimate 
disagrees. The position of those stars in the de-reddened 
color-magnitude plane suggests that they indeed have 
[M/H] > —1. In other words, for these cases [M/H]xp 
appears to be the more robust estimate. The comparison 
with GSP-Spec in the right panel of Figure 8 confirms 
the quality of GSP-Spec’s [M/H] estimates. But it also 
shows that the GSP-Spec values do not cover low [M/H] 
(by design, see Recio-Blanco et al. 2022). 

Because a star’s effective temperature (Te) can be 
very helpful in sample selection and characterization, 
we also estimate Teg using the same XP features and 
XGBoost trained on SDSS DR17 APOGEE Tog, as il- 
lustrated in Fig. 9. The XGBoost predictions based on 
XP spectra and WISE photometry match the APOGEE 
values within a median AT gg = 32 K, and with a mean 
difference of only 3 K. Figure 9 compares our [M/H] es- 
timates to APOGEE as a function of Gaia Gpp magni- 
tude, showing the slight degradation towards the faintest 
sources in the sample. 

We also compared our [M/H] estimates to the exten- 
sive set of photometric estimates for metal-poor stars 
from Chiti et al. (2021) based on SkyMapper DR2 (SM2) 
photometry. We find a mean [M/H] offset of 0.28 dex, 
and for stars with [M/H]s,;. ~ —1.5 a central 68% in- 
terval of p( [M/H]gyr2 — [M/H] xp ) of 1 dex. Clearly, the 
cross-validation of our XP-based [M/H] estimates and 
SkyMapper DR2 shows far more scatter than our cross- 
validation with external spectroscopic data sets. This 
illustrates that for this sample XP-based [M/H] _ esti- 
mates outperform those of SM2. 
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Figure 9. Top: Validation of our Te(XP) estimates against 
the Teg from SDSS DR17 APOGEE. The XP predictions 
match the APOGEE values within a median AT ess = 32 K, 
and with a mean difference of only 3 K. Bottom: Validation 
of our [M/H] estimates as a function of Ggp magnitude, with 
A[M/H] = [M/H]xp — [M/H] ,poapr, Showing that the qual- 
ity of the estimate depends somewhat on Gpp, as expected, 
but is overall unbiased. 
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